%b00(x)u0''+h00(x)u0'+b10(x)u1''+h10(x)u1'=f1
%b01(x)u0''+h01(x)u0'+b11(x)u1''+h11(x)u1'=f2
%with Dirichlet boundary 
clear all

%construct mesh
n=4000;
unimesh=(0:1/n:1)';

%construct coefficient function handel  
thebij;
load funcbij;
thehij;
load funchij;
%construct coefficient matrix A
%A=[matelt(b00,h00,n),matelt(b10,h10,n);matelt(b01,h01,n),matelt(b11,h11,n)];
A=[matelt(bij{1}{1},hij{1}{1},n), matelt(bij{2}{1},hij{2}{1},n); ...
    matelt(bij{1}{2},hij{1}{2},n), matelt(bij{2}{2},hij{2}{2},n)];

%construct source f
b=srcf(n);

%solve the equation
u=A\b;

%measure the error
uex=[u0ex(unimesh(2:n));u1ex(unimesh(2:n))];
err=max(abs(uex-u))

